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ABSTRACT 

A Fourier-Chebyshev spectral method for the incompressible Navier-Stokes 
equations is described. It is applicable to a variety of problems including 
some with fluid properties which vary strongly both in the normal direction 
and in time. In this fully spectral algorithm, a preconditioned iterative 
technique is used for solving the implicit equations arising from semi- 
implicit treatment of pressure, mean advection and vertical diffusion terms. 
The algorithm is tested by applying it to hydrodynamic stability problems in 
channel flow and in external boundary layers with both constant and variable 
viscosity. 
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INTRODUCTION 


Fourier-Chebyshev spectral methods have been employed in a number of 
numerical simulations of stability and transition in three-dimensional wall- 
bounded shear flows. Specific algorithms have been developed for straight 
channels ([1], [2], [3]), curved channels [4], the parallel boundary layer 
[5], cylindrical Couette flow [6] and pipe flow [5]. In all of these methods 
Chebyshev expansions are employed in the direction normal to the walls and 
Fourier methods are used in the remaining two directions. Hence these methods 
are applicable whenever periodic boundary conditions are appropriate in two 
directions. 

These methods usually handle the pressure and vertical diffusion terms 
implicitly, the pressure term so that the incompressibility condition is 
enforced and the vertical diffusion term in order to relax the diffusive time- 
step limitation. (The only exception is the method [4] which eliminates the 
pressure by a clever choice of divergence-free velocity expansion functions.) 
Algorithms which employ time-splitting ([1], [5]) can achieve a relaxation of 
the advective time-step limit by a semi-implicit treatment of the streamwise 
advection. These implicit equations are solved by a direct method for which 
the efficiency depends upon simple mean velocity profiles and constant 
viscosity. However, there are situations in which time-splitting errors are a 
serious problem [6]. 

If the spectral discretization in the normal direction is replaced with a 
finite difference method, then the direct solution of the Implicit equations 
can be performed efficiently for mean flow profiles and viscosities with an 
arbitrary dependence upon both the normal coordinate and time. Such Fourier- 
finite difference codes have been utilized both for channel flow [7] and for 
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the parallel boundary layer [8]. The price for this extra flexibility, 
however, is greatly reduced accuracy in the normal direction. 

The contribution of the present paper is the description of a Fourier- 
Chebyshev algorithm for wall-bounded shear flows which combines the accuracy 
and efficiency of a fully spectral scheme with the flexibility of a Fourier- 
finite difference method. The key feature of this algorithm is a precondi- 
tioned iterative technique for solving the implicit equations arising from the 
semi -implicit treatment of the pressure, mean flow and vertical diffusion 
terms. This algorithm is applicable to most of the cases described above — 
channel flow, parallel boundary layers, curved channel flow and cylindrical 
Couette flow. Relatively minor modifications are required to treat the 
different cases. Illustrations will be provided here for straight channel 
flow and for parallel boundary layer flow with constant and variable 
viscosity. The discussion will be restricted to two-dimensional flow. The 
addition of a second periodic direction is straightforward. 


2. DISCRETIZED NAVIER-STOKES EQUATIONS FOR CHANNEL FLOW 

The rotation form of the two-dimensional incompressible Navier-Stokes 
equations is 

U t ’ V(v x " V + P x = ^ U X J X + <VU y ) y (2.1) 

v t + u(v x - u y ) + p y = (yv x ) x + (y v y ) y (2.2) 


(2.3) 
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where the variable P denotes the total pressure and subscripts denote 

partial derivatives. The viscosity p is presumed to depend upon y and t 
only and the density is taken to be unity. Periodic boundary conditions in 
x and no-slip boundary conditions at y = ±1 are assumed. 

The spatial discretization of Equations (2.1) - (2.3) employs spectral 
collocation. The collocation points are 

x = j L /K j = 0, !.,••• ,K-1 (2.4) 

J X 

y m = cos (~^) m = 0,1, (2.5) 

where is the periodicity length in the streamwise direction, and K and 

N are the number of intervals in the x and y directions, respectively. 
The dependent variables have Fourier-Chebyshev series of the form 

K/2-1 N 2iTikx/L 

u(x,y,t) = l l u, (t)e X T (y), (2.6) 

k=-K/2 n=0 Xn n 

where T n is the Chebyshev polynomial of degree n. In the spectral colloca- 
tion method, spatial derivatives of u are obtained by differentiating the 
series expansion with the expansion coefficients u^ n (t) determined by 
discrete Fourier and Chebyshev transforms of the grid point values of u. The 
details of this procedure can be found in [9] and [10]. In the temporal 

discretization, the pressure gradient term and the incompressibility 

constraint are best handled implicitly. So, too, are the the vertical 

diffusion terms because of the fine mesh-spacing near the wall. The variable 

viscosity prevents the standard Poisson equation for the pressure from 
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decoupling from the velocities in the diffusion term. A simple time 
discretization uses Crank-Nicholson on the implicit terms and second-order 
Adams-Bashforth on the remainder. After a discrete Fourier transform in x, 
the following set of ordinary differential equations and boundary conditions 
result 


0 A n+l , A n+1 .r A n+l 
-gu + u + ikQ 

yy 


A n+1 . A n+1 
-gv + v 

yy 



n . A t ( , u n u n 1 \ 
= u + — l SHj " Hj J 


n-h 


« A A 

- ikQ + gu 

yy 



(2.7) 

( 2 . 8 ) 


and 


AA n+l A n+1 n 
-iku - v = 0, 

y 


u(-l) = u(+l) 


(2.9) 


( 2 . 10 ) 

A A 


v(-l) - v(+l) - 0 

A . | 1 *. A * f- * 

where k = 2irk/L x> 3 = + > Q “ ~ P, i = /-l, and hats denote Fourier 

transformed variables in wavenumber space* The wavenumber is denoted by k 

A A 

and the dependence of u, v, and Q upon k has been suppressed. The 
superscript n represents the time level. and H 2 , which contain the 

terms treated explicitly, are given by 


H, = -v(u - v ) + (pu ) + p u - P I 

1 y x x x y y x'mean 


( 2 . 11 ) 


H„ = -u(v - u ) + (pv ) + p v . 

2 x y xx y y 


( 2 . 12 ) 


The last term in Eq. (2.11) is the mean streamwise pressure gradient which 
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drives the channel flow. All of these derivatives are evaluated by spectral 
collocation. A semi-implicit treatment of the mean streamwise advection term 
is easily incorporated. For example, the left-hand side of Eq. (2.7) has the 
additive term 



0 


"n+1 
u ; 


in addition, ug u x appears in Eq. (2.11). 

For each wavenumber k, the system of Eqs. (2.7) - (2.9) can be written as 


L U » F, (2.13) 

where U = (u 11 *^, v 1 *^, Q n+ ^) and F is the known right-hand side. The 
matrix L is a full MxM matrix where M = 3N. A direct solution of (2.13) 
by Gauss elimination methods would require 0(M Z ) storage and 0(M J ) 
arithmetic operations. An iterative solution, on the other hand requires only 
0(M) storage and 0(M log M) operations per iteration. The description of 
an effective iterative scheme will be provided in the next section. The use 

Ai A 

of the variable Q in place of P puts L into a nearly self-adjoint form. 

At this point some remarks pertinent to our selection of this scheme are 
in order. Our goal was to develop a single, fully spectral algorithm which is 
applicable to a broad class of problems. Many interesting phenomena involve a 
strong variation of the viscosity, the mean advection, and/or the geometric 
terms in the direction normal to the wall (or walls) and possibly also in 
time. (A number of three-dimensional calculations employing the present 
algorithm on such problems are in progress and will be reported elsewhere.) 
In many of these problems semi-implicit treatment of the normal diffusion 
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and/or the mean strearawise advection are desirable. The observations of 
Marcus [6] about the pitfalls of time-splitting in some problems is a strong 
argument in favor of an un-split method for a general purpose algorithm. A 
Chebyshev tau method in the normal direction is ruled out in favor of 
Chebyshev collocation in all but the simplest cases. The variable viscosity 
and mean advection prevent the velocity and pressure equations from de- 
coupling as in the influence matrix methods ([3], ([6]). The matrix 
diagonalization technique for solving Eq* (2*13) is not practical because the 
matrix L may depend upon time. These considerations have led us to develop 
an iterative technique for solving the collocation equations. 


3. SPECTRAL SOLUTION WITH FINITE DIFFERENCE PRECONDITIONING 

The key to the efficiency of an iterative method for the solution of Eq. 
(2.13) is the use of an effective preconditioning matrix so that the number of 
iterations is small. The reason is that the condition number of the matrix 
L is large. Consequently, standard iterative techniques would be slow. But 
let H be some preconditioning matrix for L, i.e. , the iterative scheme is, 
in effect, applied to the equation 

H -1 L U = H' 1 F. 


The desirable properties of the preconditioning matrix are that the condition 
number of H“* L be small and that equations such as 


H U = G 
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can be solved cheaply for U (relative to the evaluation of LU). The first 
property implies that only a small number of iterations are required and the 
second property implies that a single preconditioned iteration costs roughly 
the same as a single un-preconditioned iteration. We base our choice of H 
on Orszag's suggestion [11] that a finite difference approximation to the 
differential equation be used. The interesting physical problems have high 
Reynolds number, i.e., low viscosity. Thus, the first derivative terms in Eqs. 
(2-7) - (2-9) predominate. Therefore, the effective preconditioning of them 
is crucial. 

To illustrate the difficulty with first derivative terms and to assess 
various remedies we consider the model scalar problem 

u^ = f (3.1) 

on [0,2ir] with periodic boundary conditions. The appropriate spectral 
method uses Fourier collocation. The eigenfunctions of the discrete spectral 
operator L and of the finite difference operator H are the exponentials 

ikx. 

« 3 


where k is the wavenumber and Xj is a Fourier collocation point as given 
by Eq. (2.4). Four possibilities for the finite difference operator are 
considered here: central differences, central differences with a high mode 

cutoff, one-sided differences and the use of a staggered mesh. The 

eigenvalues of these preconditioned matrices, H * L, for the model scalar 
problem are given in Table I for all four possibilities. The term kAx is 
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the product of the wavenumber k and the grid spacing Ax. It falls in the 
range 0 < kAx < 2ir . The eigenvalues for the centered differences 
kAx/sin kAx, are unbounded as kAx > tt. Thus, pure central difference 
preconditioning yields a large condition number for L. Orszag [11] noted 

that truncating the high modes limits the eigenvalues. Table I indicates that 
this does produce a bounded spectrum; the price is that some high wavenumber 
information is lost. Another cure is to use one-sided (forward or backward) 
differences for the first derivative terms. For the model problem, it results 
in bounded but complex eigenvalues with real parts tending to zero. Many 
iterative schemes perform badly on such problems. For the staggered mesh the 
eigenvalues of the preconditioned matrix for the model problem remain bounded 
and real, with no loss of high wavenumber information. 

These model problem results led us to consider a staggered mesh for the 
Navier-Stokes equations. The staggered mesh which is appropriate for the 
Fourier^Chebyshev discretization is staggered only in the y direction. The 
velocities are defined at the cell faces y m , as given by Eq. (2.5), and the 
pressure is defined at the cell centers 

y m _ y = cos(ir(m - V 2 )/ N ) m = (3.2) 

The momentun equations are enforced at the faces, whereas the continuity 
equations are enforced at the centers. The velocity boundary conditions are 
enforced at the two walls. Note that there is no need for an artificial 
pressure boundary condition at the walls. 

The staggered mesh assigns one less vertical degree of freedom to the 
pressure than to the velocities. This is common practice in finite element 
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techniques for the Navier-Stokes equations (see, for example, [12]). Huberson 
and Morchoisne [13] have recently proposed a filtering procedure for spectral 
solutions of the incompressible Navier-Stokes equations on a non-staggered 
mesh. It has the effect of removing one vertical degree of freedom from the 
pressure. 

Let us now examine some of the mechanics involved in employing a staggered 
mesh for Eqs. (2.7) - (2.9). Focus first on the spectral evaluation of the 
various terms. The explicit terms, denoted by and H 2 , are evaluated in 
a straightforward manner since they are required at the faces and involve only 
the velocities. The same holds for the remaining velocity terms in the 
momentum equations. The only complication here is the two terms involving the 
pressure. From the values of Q at the centers, trigonometric interpolation 
can be used to obtain Q at the faces. First, use the center values to 
obtain the Chebyshev coefficients 


P 

n 


i I *W 2 ) 


cos 


irn(m - Vg ) 


N 


m 


(3.3) 


for n = 0,1, ### ,N-1, where the dependence upon k and t has been 

A A 

suppressed. Then set = 0 and compute the values of Q at the faces 


N , 

Q(y m ) = E On COS 

n=0 


nmn 

N 


m = 0, 1 ,• • • ,N. 


(3.4) 


Both of these sums may be computed by fast cosine transforms. This takes care 
of the pressure term in Eq. (2.7). The term in Eq. (2.8) may be 

A 

evaluated from the values of Q at the faces in a standard fashion. For the 
continuity equation one first evaluates 


r = ik u + Vy 

at the cell faces In the standard manner and then interpolates this result to 
the cell centers, via 


r 

n 



Nc eofO 
n 



r(y ) cos 
w m 


irnm 

N 


for n * 0,1,»*»,N where 


2 n = 0 or N 

9 

1 1 < n < N-l 


(3.5) 


(3.6) 


and 


* N 

r( W Jo 


r cos 
n 


(m~ V 9 )n 
N 


m 


1 >* * * ,N. 


(3.7) 


The finite difference operator H pertains only to the left-hand side of 
Eqs. (2.7) - (2.9). The second derivative of the velocities is evaluated by 
3-point centered differences of the values at the faces, using the formula 
appropriate for the non-uniform grid, e.g. , 


yy 


m 


2 r m+1 m 

^nri- 1 y m- 1 y mt 1 y m 


U m ~ Vli 

y m “ Vl 


The pressure term in the u momentum equation is approximated by a linear 
average of the adjacent cell— centered pressure values. The vertical pressure 
gradient term in the v momentum equation is approximated by 2 — point 


differences of the adjacent cell— centered pressure values. The streamwise 
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velocity in the continuity equation is taken as the linear average of the 

A 

velocity values at adjacent cell faces and the vertical derviative of v uses 
2-point differences of the adjacent cell faces values. 

Order the unknowns as 

u ■ ( u !» v !> Qi/2’ u 2’ V 2’ ^3/2 ** * * * U N’ V N’ ^N-l/2^ 
and order the equations as 

continuity at y^/2 
v momentum at 
u momentum at y^ 


continuity at Yn- 3/2 
v momentum at 
u momentum at y^-l 

u BC at y N 
v BC at y N 

continuity at yN-1/2* 

This requires the velocity boundary conditions at yQ to be absorbed into the 
matrix. This ordering produces a block tridiagonal structure for H that can 
evidently be solved without pivoting within the diagonal block. (We have no 
proof for this claim, but we have made numerous checks. In all cases the 
solution without pivoting produced results that agreed with solutions with 


pivoting to at least 8 digits.) 
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For k = 0 the structure is even simpler. The velocity v is constant 
in y, the velocity u satisfies a tridiagonal equation, and the pressure 

A 

Q satisfies a bidiagonal equation. The latter is solved by setting 


Q(y ij ) = 0 and then solving for each successive value of pressure. This 

A 

particular choice of Q(y y ) is arbitrary and corresponds to specifying the 
mean level of pressure. 

We have computed eigenvalues of L not only for the staggered grid 


method but also for the same three alternatives that were discussed for the 


model problem. In these cases the pressure is defined at the cell faces and 
the continuity equation is enforced at the cell faces. This version requires 
numerical boundary conditions for the pressure at the walls. The continuity 
equation and the vertical momentum equation yield 


P = -ik(yu) . 

y y 

The finite difference approximation uses one-sided differences and the 
matrix H is still block tridiagonal. 

The eigenvalues of the preconditioned matrix H * L for the Navier-Stokes 
equations are displayed in Figures 1 and 2 for two wavenumbers and for four 
different discretizations of the first derivative terms. The calculations 

were made for a K = 32, N = 16 grid, with y = (7500)“* and a streamwise 
CFL number of 0.1. 

A 

The results for k = 1 are particularly interesting. When central 
differences for the first derivative terms are used, there are several complex 
eigenvalues with large real parts. The remaining eigenvalues are real with 
1.0 < A < 4.5. As N increases, both the real and imaginary parts of the 
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eigenvalues grow. (The largest eigenvalues for N = 24 and 32 are 12.3 ±i 
4.5 and 16.5 ±i 6.4 respectively.) When the upper one-third of the Chebyshev 
modes are cut off in the first derivative representation the spectrum is 
apparently bounded from above. However, there are now a number of complex 
eigenvalues with small real parts. One-sided first differences yield mainly 
complex eigenvalues including some with very small positive real parts. When 

A 

the mesh is staggered, all the eigenvalues for k = 1 lie close to the real 
axis between 1 and ir/2 = 1.57, 

The eigenvalue spectra are only slightly different at higher wavenumbers, 

A 

as illustrated in Figure 2 for k = 10. Although there are some complex 
eigenvalues for the staggered mesh, they are reasonably well-confined. 
Similar eigenvalue calculations have been performed for the staggered grid 
algorithm for N = 24 and N = 32. The real parts are still confined 
between 1 and tt/2 and the magnitudes of the imaginary parts decrease as 
N increases. 

Note that the model problem estimates the eigenvalue trends surprisingly 
well considering that it is just a scalar equation, has only first derivative 
terms, and uses Fourier series rather than Chebyshev polynomials. 

The preceding results indicate that the staggered grid leads to the most 
effective treatment of the first derivative terms. The condition number of 
the preconditioned system is reasonably small and full resolution is retained. 
However, the iterative scheme used for solving Eqs. (2.7) - (2.9) must be 
capable of dealing with the complex eigenvalues. Two types of iterative 
schemes are feasible. Chebyshev iteration [14] will converge because the real 
parts of the eigenvalues are greater than 1. However, this method contains 
parameters that depend upon the location of the eigenvalues in the complex 



14 


plane. Alternatively, a parameter-free variational method [15] such as the 
minimum residual (MR) method will work provided that the Hermitian part of 
LH - 1 is positive definite. This condition is satisfied for all the cases 
discussed in this paper. 

The preconditioned version of MR for Eq. (2.13) involves making an initial 
guess U®, computing the initial residual 


R° = F - LU° , 


solving 


HZ° = R° 


and then iterating according to 


= (LZ £ , R £ ) 

1 (lz £ ,lz £ ) 


U £+1 = U £ + a Z £ 
£ 


R £+i = r £ L Z £ 


Hz £+1 - R £+1 


(3.8) 

(3.9) 

(3.10) 

(3.11) 

(3.12) 

(3.13) 


until convergence. The parameter in Eq. (3.10) is chosen so that the 

residual in Eq. (3.12) is as small as possible consistent with the 
prescription (3.11). Representative convergence histories for the MR method 
are shown in Figure 3 where the Lj norm of the residual is plotted against a 
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number of iterations for N = 16, 32, and 64. At a Reynolds number of 7500, 
each iteration is found to reduce the residual by almost an order of magnitude 
and there is a trend of faster convergence with increasing N which may be 
partly attributed to the higher resolution. The physical results to be 
presented in Figure 5 and Table III become insensitive when the Lj norm is 
smaller than 10“^. 

4. EVOLUTION OF SMALL DISTURBANCES IN CHANNEL FLOW 

In order to test the algorithm proposed for Navier-Stokes equations, we 
study the problem of the evolution of small disturbances in channel flow. 
This problem has been studied extensively using the Orr-Sommerfeld equation. 
When the amplitude of the disturbances imposed upon the mean (time 

O 

independent) channel flow u(y) = (1 - y c ) is small, then the numerical 
solution of the Navier-Stokes equation should be the same as that implied by 
the Orr-Sommerfeld solution. This linear solution has the form 

u(x,y,t) = (1 - y 2 ) + e Re{ip y (y) , (4.1) 

v(x,y , t) * -e Re{ia\Ky)e i ^ ax-(1Jt ^} , (4.2) 

where is the eigenfunction normalized to a maximum value of 1, tu is the 
complex frequency (with the largest imaginary part), a is the prescribed 
wavenumber, and e is the perturbation amplitude. 

The perturbation flow energy E(t) is 

L x 1 

E(t) = / d x / { [u(x,y,t) - ( 1— y 2 ) ] 2 + v 2 (x,y,t)}dy, 

0 -1 


(4.3) 
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where L = 2tt/<x. Choose initial conditions from Eqs. (4.1) and (4.2) with 

2o) . t 

t = 0 and let Eq = E(0). For small amplitudes E(t)/EQ = e 

The particular problem chosen for study had y = (7500)“* and a ■ 1. The 
only unstable mode has a) = 0.24989154 + i. 00223498. The amplitude parameter 
was e = .0001. Two different discretizations in y were used: (1) 
Chebyshev collocation and (2) finite differences. Both methods used a Fourier 
collocation method in x. The Fourier-finite difference method used a 
staggered mesh, with the cell centers given by Eq. (3.2) and the cell faces 
located midway between the neighboring cell centers. This method is just that 
of Moin and Kim [7], applied to a direct simulation. Only four collocation 
points were used in the x-direction. For this basically linear test problem, 
the x-direction has essentially perfect resolution. The time step was small 
enough so that the vertical discretization errors were predominant in all but 
the most highly resolved cases. 

The basic comparison of the vertical discretizations is provided in 
Figures 4 and 5, where the natural logarithm of the perturbation energy ratio 
is plotted. The solid line in the figures represents the linear stability 
result. The finite difference solution is plotted in Figure 4 for several 
vertical grids. Even the N = 256 results are appreciably in error. 

The Fourier-Chebyshev results are presented in Figure 5. The results for 
the N - 32 grid are already in excellent agreement with the linear theory 
results. The numerical results for N = 16 are wildly inaccurate. This is in 
contrast with the finite difference calculations where £n E/Eq at least 
varies linearly with time for various grids. This behavior is typical of 
spectral methods in general: if the resolution is inadequate, say worse than 
20%, then the spectral results are inferior to finite difference results; 
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however, once the 10% accuracy level is achieved, spectral results become 
dramatically superior. 

In the above calculations, all runs were terminated at t/tg = 2 where 
tg is the time required for the wave to propagate through the streamwise 
computational domain. In this case, tg = 25.1438. The calculated energy 
ratio and its error at one and two periods are presented in Tables II and III 
for the finite difference and the Chebyshev methods, respectively. The 
convergence of the finite difference method is quadratic. The convergence of 
the Chebyshev method is dramatic: the N = 32 spectral results are far better 
than N = 256 finite difference results (and took less CPU time). The error 
for the N * 64 Chebyshev case is dominated by time discretization and 
nonlinear effects. 

The spectral results were all obtained with a time-step corresponding to a 
mean streamwise CFL number of 0.025 and with an explicit treatment of 
advection. Such a small time-step is necessary for accuracy purposes. 
Stability problems over two periods only arise for CFL numbers above 0.30. 
The advantage of the capability of the algorithm to treat the mean advection 
implicitly arises in calculations with higher spatial resolution. An example 
is provided by calculations for this same test problem using 16 Fourier modes 
rather than 4. The semi-implicit advection version of the algorithm is stable 
for CFL numbers as large as 1. However, the accuracy suffers for such large 
time-steps. 
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5. DISCRETIZATION FOR EXTERNAL BOUNDARY LAYERS (0 < y < r^) 


This numerical method may also be applied to a model of the external 
boundary layer. In order to use periodic boundary conditions in the 

streamwise direction, one must make the parallel flow assumption, i.e. , fix on 
some reference location in a spatially growing boundary layer and use the 
corresponding mean velocity profile at all x. One must then set the mean 
vertical velocity to zero and make a minor adjustment to the mean streamwise 
pressure gradient to achieve parallel flow. 

A stretching transformation can be applied in the (unbounded) vertical 
direction. Let 


y 


a 


1 + S 
b - K 


(5.1) 


where y is the physical vertical coordinate, £ is the computational 
coordinate and a and b are constants. Let be the upper boundary in 
the physical plane and set 

b = 1 + — . (5.2) 

n„ 


Then for any choice of the scaling parameter a, the computational coordinate 
5 falls within the standard Chebyshev interval [-1,1]. Derivatives in the 
vertical direction are evaluated by multiplying the Chebyshev collocation 
derivative in £ by the Jacobian of the transformation, i.e., 



(5.3) 


The necessary modifications to Eqs. (2.7) to (2.12) are straightforward 
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A number of choices are available for the numerical boundary condition at 
T^. The simplest is to require that the solution at y - correspond to 

the flow at infinity. This is accomplished by setting u at y = for 

k = 0 equal to the free stream velocity and setting all other velocity 
components to zero. Another approximation was used by Fasel [16] in his 
finite difference calculations of the boundary layer: 



(5.4) 

A A A 

v y = - |k|v. 

These two alternatives will be referred to below as the zeroth-order and 
first-order boundary conditions, respectively. 

The finite difference preconditioning matrix is straightforward. Both 
types of upper boundary conditions lead to a block tridiagonal structure for 
H (which does not appear to require pivoting). The eigenvalues for the pre- 
conditioned matrix are illustrated in Figure 6. The grid has K * 32, N = 16, 
and the Reynolds number (y * ) is 7500. Three different CFL numbers (0.01, 
0.1, 1.0) are checked and the effect is found to be negligible. The eigen- 
values do tend to become widely apart with increasing ri^. For fast 
convergence, therefore, one would like to impose the freestream boundary 
conditions at as small q as possible. Representative convergence histories 
of the MR method for the boundary layer case are shown in Figure 7. Boundary 
conditions are imposed at = 10 and 20. Both the zeroth and first-order 
boundary conditions are used and the convergence is found to be significantly 
faster in the latter case. The physical results to be presented become 
insenstive when the Lj norm is smaller than 10"^. 
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We now describe results of computations of the evolution of small 
disturbances in flat plate flow with no-slip boundary conditions at the solid 
wall. The initial conditions are the Orr-Sommerfeld solution imposed upon the 
Blasius profile. 

The particular problem chosen for study had y = (1500) - * , a = .3 and 
a) = .10288548 + i. 00249003. The zeroth-order boundary conditions were 

imposed at ^ = 20. The amplitude parameter e is taken to be .0001. Four 
Fourier spectral modes were used in the streamwise direction. All runs were 
terminated at t/tg = 2 where tg is the time required for the wave to 
propagate through the streamwise grid. In this case tg = 61.0689. 

Results analogous to those provided earlier for channel flow are given in 
Figures 8 and 9 and Tables IV and V. The N = 32 Chebyshev results are far 
more accurate than the N = 256 finite difference results and required less 
CPU time. Some additional Fourier-Chebyshev calculations were performed to 
assess the upper boundary conditions. The results are reported in Table VI in 
terms of the energy error after two periods. For = 10, first-order 

boundary conditions provide a significant inprovement in accuracy over the 
zeroth-order ones. At = 20, however, the improvement is marginal. The 

results for first-order boundary conditions at n = 10 are plotted in Figure 
10. Significant improvement for N = 16 can be noted in comparison with 
Figure 9 where zeroth-order boundary conditions were imposed at n a = 20, 

In order to test the variable viscosity capability of the numerical 
algorithm, we applied it to water boundary layers with wall heat transfer. 
The viscosity of water is a strong function of temperature, decreasing with 
increasing temperature. Thus, heating of water boundary layers has a 

stability effect. We used the empirical temperature— viscosity formula given 
in [17] . 
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An Orr-Sommerfeld equation for incompressible flow to include the effect 
of viscosity can be derived as in [18] by neglecting temperature perturba- 
tions. This equation has been solved to provide initial conditions for the 
Navier-Stokes code which also neglects the temperature perturbation. The free 
stream temperature 0^ is assumed to be 293° K and three wall-to-f reestream 
temperature ratios were examined: ® w /®oo = 1»0, 0.9. The resulting 

viscosity distributions calculated are plotted in Figure 11. The f reestream 
Reynolds number (p^) ^ = 10000 and a = 0.15. The Orr-Sommerfeld 

eigenvalues and time periods (tg) for the three cases are given in Table 

VII. Hie Navier-Stokes solutions for the three temperature ratios are 

presented in Figure 12. The solution has been obtained with first-order 
boundary conditions imposed at q = 20 and using N = 32. In each case 500 
time steps were used to reach t/tQ = 2. The solid line in each case 

represents the theoretical results. While the growth rates are vastly 
different for the three cases, the calculated results are in good agreement 
with the theory. The calculated energy ratios and errors are given in Table 

VIII. 

In Figure 12, a strong stabilization effect may be noted when the wall-to- 
f reestream temperature ratio 0^/0^ is increased from 1 to 1.1. These 
calculations were performed with £ = .0001 using four Fourier spectral modes 
in the streamwise direction. In order to study the effect of nonlinearity, we 
have recomputed the = case using eight Fourier spectral modes in 

the streamwise direction for e = .0001, .01, .03, .05. The results are 

presented in Figure 13 along with the linear Orr-Sommerfeld solution. It can 
be seen that the energy rate increases with increasing perturbation amplitude 
e. A thorough set of two-dimensional and three-dimensional finite amplitude 
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results, produced by the latter two authors in collaboration with D. Bushnell, 
will be presented elsewhere. 


6. CONCLUDING REMARKS 

A Fourier-Chebyshev spectral method for the solution of the incompressible 
Navier-Stok.es equations has been presented. This fully spectral method is 
applicable to both the internal and external boundary layers with variable 
viscosity. The method uses Chebyshev polynomials in the vertical direction 
and Fourier spectral collocation in the horizontal direction. The continuity 
and momentum equations are solved as a set of coupled equations without 
splitting. A staggered grid is employed in the vertical direction so that no 
numerical pressure boundary conditions are needed. The resulting implicit 
equations are solved by a preconditioned iterative technique. The algorithm 
has been subjected to extensive testing by applying it to problems in 
hydrodynamic stability in channel flow and external boundary layers with 
constant and variable viscosity. The results obtained with 33 Chebyshev 

polynomials are found to be much more accurate and require less CPU time than 
when 257 finite difference grid points are used. 
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Table II. Channel Fourier-Finite Difference Convergence 


N 

1 period 
Ef/^ol calc* 

E f/ E 0 1 error 

16 

.31369085 

-.80526006 

32 

.59348926 

-.52546165 

64 

.93539768 

-.18355323 

128 

1.06837752 

-.05057339 

256 

1.10598936 

-.01296155 


N 

2 periods 
Ef /^ol calc. 

E f/ E oUrror 

16 

.45883275 

-.79321839 

32 

.26725477 

-.98479637 

64 

.81641093 

-.43564021 

128 

1.12221673 

-.12983441 

256 

1.21820807 

-.03384307 
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Table III. Channel Fourier - Chebyshev Convergence 



1 period 


N 

E f 1 calc. 

E f/ E 0 1 error 

16 

1.17188803 

.05293712 

32 

1.11912239 

.00017148 

64 

1.11896735 

.00001644 


N 

2 periods 
E f ^ E 0 1 calc. 

E f/ E 0 lerror 

16 

2.07329163 

0.82124050 

32 

1.25291992 

.00086879 

64 

1.25214542 

.00009429 
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Table IV. Boundary Layer Fourier-Finite Difference Convergence 


N 

1 period 
E^/Eq ( calc. 

Ef / E oUrror 

32 

1.4144405 

0.05899918 

64 

1.2294246 

-0.12601673 

128 

1.3213179 

-0.03412346 

256 

1.34699412 

-0.00844721 


N 

2 periods 

V E olcalc. 

E f /E ol error 

32 

4.4395339 

2.6023126 

64 

1.5376562 

-0.29956505 

128 

1.7459684 

-0.09125289 

256 

1.8144658 

-0.02275544 
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Table VI. Effect of ri^ and Top Boundary Condition 


^f/^olerror 


N 


= 10 

oo 

= 20 


0th order 

1st order 

0th order 

1st order 

16 

-0.21312422 

0.16289710 

-0.51718473 

-0.52723330 

32 

-0.09173191 

0.00023775 

0.00048000 

0.00021268 


Table VII. Orr-Sonmerfeld Solution for Hater Boundary Layer 
with Hall Heat Transfer (a = 0.15, (p w ) * = 10000) 


6 /e 

w <® 

0) 


t 0 

1.1 

0.02872049 + i 

0.00020520 

218.7701 

1.0 

0.03386607 + i 

0.00343206 

185.5303 

0.9 

0.03445962 + i 

0.01259238 

182.3347 
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Table VIII. Navier-Stokes Solution for Water Boundary Layer with 

Wall Heat Transfer (y = 20, N = 32, 1st order B.C.) 

'max * 



1 period 


0 /0 
W <» 

E f /E olcalc. 

E f /E 0 1 error 

1.1 

1.0946078 

0.00067076 

1.0 

3.5720516 

-0.00127094 

0.9 

99.374109 

0.67521093 


0 /0 
W «> 

2 periods 
E f /E olcalc. 

E f ^ E 0 1 error 

1.1 

1.1972253 

0.00052696 

1.0 

12.757289 

-0.01134543 

0.9 

9871.3827 

129.91024 
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Figure 1. 


Figure 2. 


Figure 3. 


Figure 4. 


Figure 5. 


Figure Captions 

A plot of k = 1 channel flow eigenvalues of the preconditioned 
matrix H - ^ L for four first derivative treatments. In this case 
y -1 = 7500, CFL = 0.1, N = 16, and K = 32. 

a 

A plot of k = 10 channel flow eigenvalues of the preconditioned 
matrix L for four first derivative treatments In this case 

p -1 = 7500, CFL = 0.1, N = 16, and K = 32. 

Convergence history of the minimum residual method for the channel 
flow problem (p 1 ■ 7500), 

Computed perturbation energy ratio for channel flow problem 
(p ^ = 7500). A Fourier spectral method in x and a second-order 
finite difference method in y are used. Results are shown for a 
four point grid in x and for various grids in y. The solid 
line is the correct result. 

Computed perturbation energy ratio for channel flow problem 
(v ^ ** 7500). A Fourier spectral method in x and a Chebyshev 
spectral method in y are used. Results are shown for a four 
point grid in x and for various grids in y. The solid line is 


the correct result 
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Figure 


Figure 


Figure 


Figure 


Figure 


6. A plot of the eigenvalues of the preconditioned matrix H~^ L for 
an external boundary layer using the staggered grid. In this 

case, y * = 1500, N = 16, and K = 32. Zeroth-order boundary 
conditions are imposed at n . 

7. Convergence history of the minimum residual method for the 
boundary layer problem. Zeroth-order boundary conditions are used 
in the parts of the figure on the left-hand side and first-order 
conditions in parts on the right-hand side. 

8. Computed perturbation energy ratio for the boundary layer problem 

with constant viscosity (y * = 1500). A Fourier spectral method 
in x and a second-order finite difference method in y are 
used. Results are shown for a four-point grid in x and various 

grids in y. The solid line is the correct result. 

9. Computed perturbation energy ratio for the boundary layer problem. 

A Fourier spectral method in x and a Chebyshev spectral method 
in y are used. Results are shown for a four-point grid in x 

and various grids in y. The solid line is the correct result. In 
this case, zeroth-order boundary conditions are imposed at 

n = 20. 

oo 

10. Computed perturbation energy ratio for the boundary layer problem. 

A Fourier spectral method in x and a Chebyshev spectral method 
in y are used. Results are shown for a four-point grid in x 
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and two different grids in y. The solid line is the correct 

result. In this case, first-order boundary conditions are imposed 

at n = 10. 

00 

Figure 11. Variation of viscosity for a water boundary layer with and without 

wall heat transfer (0 = 293° K). 

00 

Figure 12* Computed perturbation energy ratio for a water boundary layer 

= 10000) using a Fourier-Chebyshev spectral method. The 
results shown are for a four-point grid in x and a 33 point-grid 
in y. 0^/0^ = 1*1 > 1.0, 0.9 pertain to wall heating, no heating 
and wall cooling respectively. The solid lines are the correct 
results. 

Figure 13. Computed perturbation energy for a water boundary layer 

(y 10000, 0 /0 = 1.1). The results are shown for various 

initial perturbation amplitudes to indicate the effect of 
nonlinearity and were computed by using an eight-point grid in x 
and a 33-point grid in y. The solid line is the linear result. 
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